
clear all 
folderfig   = '/Users/wittwer/Dropbox/5_Jason and Milena/data/Matlab_JMP_2R1/main_estimation_submission2R1/3_figures';
foldermats  = '/Users/wittwer/Dropbox/5_Jason and Milena/data/Matlab_JMP_2R1/main_estimation_submission2R1/2_output';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%CHOOSE SPECIFICATION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%A) benchmark model -> for the model fit

    %largetradeL=0;
    %tradesizeL =0;
    %loyaltyL=0; and =1

%B) extended model with small quantities -> for the main counterfactual analysis

    %largetradeL=1;
    %tradesizeL =0;
    %loyaltyL=0; and =1
    
%For both specifications pick: 

    %zerobenefit =0;
    %homedealerL=1;
    %dayL=1;
    %keepallL=1;
    %fejL=1;
    %kappaL =0.1;
    %etaL=0; 
    %quantityrobusL=0;

%%choose whether to show plots 
plot_figuresL=0;

zerobenefit =0; %--> CAREFUL: most parameters are automatically chosen below
    %=1 only if using the homedealer model but shutting off the loyalty benefit - benchmark model with fixed base qualitites
    %=0 else
    
%%choose which model specification
homedealerL=1;
    %=0 benchmark model
    %=1 model with homedealer
      
loyaltyL =0;    
    %=0 set r=0
    %=1 keep r>0

dayL=1;
    %=0 for day (benchmark model)
    %=1 for week (extended model)
    
keepallL=1;    
    %=1 keep all invesotrs
    %=0 keep investors with LEIs and retailers
     
fejL=1; %THIS USED TO BE 1 (not exist)
          
largetradeL=0; %NOTE: for retailers always used largetradeL=0 (this is already coded below)
    %=0 all trade sizes 
    %=1 split into trade sizes (large and small)

tradesizeL =0;
    %=0 small trade
    %=1 large trade (>25 million)
    
%set to 0 if investor gets 0 rent in bilateral trade
if homedealerL==1
        xij_for_investorL=0;  
    else 
        xij_for_investorL=1;  
end
    %When homedealer=0,
        %xij_for_investor=1 means that investor earns 0 surplus in bilateral --> BENCHMAKR
        %xij_for_investor=0 means that investor earns xij in bilateral trade 
 
    %When homedealer=1,
        %xij_for_investor=1 means that investor earns 0 surplus in bilateral trade 
        %xij_for_investor=0 means that investor earns delta (home dealer benefit over quality) in bilateral trade --> BENCHMAKR

kappaL =0.1; %share of retailers
    %=0.1 for all CF but CF=7 
    %=0 for CF=7 (eliminate platform access for institutionls assuming away retail investors)
       
%%choose eta 
etaL=0; 
    %=this is the only option for this file
    
quantityrobusL=0;


if zerobenefit ==1 %for CF use the small trades only 
    homedealerL=0; dayL=1; keepallL=1; largetradeL=1; tradesizeL=0; quantityrobusL=0;
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%LOAD DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%A) IB (group 1) 
if zerobenefit==0
    temp=['/estimation',['_',num2str(1),'_' num2str(0),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(tradesizeL), '_', num2str(fejL)]]; 
    load(fullfile([foldermats,temp,'.mat']))
else 
    temp=['/estimation',['_',num2str(1),'_' num2str(0),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(tradesizeL), '_', num2str(fejL), '_zerobenefit']]; 
    load(fullfile([foldermats,temp,'.mat']))
    homedealerL=0;
end 

%attach Mhat to Y - bring in long format
Mhatlong=NaN*ones(size(Y,1),3);
for t=1:Tg
    Tt = length(new_period(t,2):new_period(t,3));
    Mhatlong(new_period(t,2):new_period(t,3),:)= ones(Tt,3).* Mhat(t,:);
end 

Y=[1*ones(size(Y,1),1),Y,Mhatlong];
vdjSQ_IB_correct=vdjSQ;

YIB=Y;
vecIB=vec;
keep YIB folderfig foldermats vdjSQ_IB_correct  vecIB xij_for_investorL plot_figuresL homedealerL etaL quantityrobusL dayL kappaL keepallL largetradeL tradesizeL zerobenefit fejL loyaltyL


%B) IS (group 3)
if zerobenefit==0
    temp=['/estimation',['_',num2str(-1),'_' num2str(0),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(tradesizeL), '_', num2str(fejL)]]; 
    load(fullfile([foldermats,temp,'.mat']))
else 
    temp=['/estimation',['_',num2str(-1),'_' num2str(0),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(tradesizeL), '_', num2str(fejL), '_zerobenefit']]; 
    load(fullfile([foldermats,temp,'.mat']))
    homedealerL=0;
end 

%attach Mhat to Y - bring in long format
Mhatlong=NaN*ones(size(Y,1),3);
for t=1:Tg
    Tt = length(new_period(t,2):new_period(t,3));
    Mhatlong(new_period(t,2):new_period(t,3),:)= ones(Tt,3).* Mhat(t,:);
end 

Y=[3*ones(size(Y,1),1),Y,Mhatlong];
vdjSQ_IS_correct=vdjSQ;
YIS=Y;
vecIS=vec;
keep YIB YIS folderfig foldermats vdjSQ_IB_correct vdjSQ_IS_correct  vecIB vecIS xij_for_investorL plot_figuresL homedealerL etaL quantityrobusL dayL kappaL keepallL largetradeL tradesizeL zerobenefit fejL loyaltyL


%C) RB (group 2) 
if zerobenefit==0
    temp=['/estimation',['_',num2str(1),'_' num2str(1),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(0), '_', num2str(fejL)]]; 
    load(fullfile([foldermats,temp,'.mat']))
else 
    temp=['/estimation',['_',num2str(1),'_' num2str(1),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(0), '_', num2str(fejL), '_zerobenefit']]; 
    load(fullfile([foldermats,temp,'.mat']))
    homedealerL=0;
end 

%attach Mhat to Y - bring in long format
Mhatlong=NaN*ones(size(Y,1),3);
for t=1:Tg
    Tt = length(new_period(t,2):new_period(t,3));
    Mhatlong(new_period(t,2):new_period(t,3),1:2)= ones(Tt,2).* Mhat(t,:);
end 

Y=[2*ones(size(Y,1),1),Y,Mhatlong];
YRB=Y;
vecRB=vec;
keep YRB YIB YIS folderfig foldermats vdjSQ_IS_correct vdjSQ_IB_correct vecIB vecIS vecRB xij_for_investorL plot_figuresL homedealerL etaL quantityrobusL dayL kappaL keepallL largetradeL tradesizeL zerobenefit fejL loyaltyL


%C) RS (group 4) 
if zerobenefit==0
    temp=['/estimation',['_',num2str(-1),'_' num2str(1),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(0), '_', num2str(fejL)]]; 
    load(fullfile([foldermats,temp,'.mat']))
else 
    temp=['/estimation',['_',num2str(-1),'_' num2str(1),'_',num2str(dayL), '_', num2str(homedealerL), '_', num2str(keepallL),'_', num2str(quantityrobusL),'_', num2str(largetradeL), '_', num2str(0),  '_', num2str(fejL), '_zerobenefit']]; 
    load(fullfile([foldermats,temp,'.mat']))
    homedealerL=0;
end 
   
%attach Mhat to Y - bring in long format
Mhatlong=NaN*ones(size(Y,1),3);
for t=1:Tg
    Tt = length(new_period(t,2):new_period(t,3));
    Mhatlong(new_period(t,2):new_period(t,3),1:2)= ones(Tt,2).* Mhat(t,:);
end 

Y=[4*ones(size(Y,1),1),Y,Mhatlong];

YRS=Y;
vecRS=vec;
keep YRS YRB YIB YIS folderfig foldermats vdjSQ_IS_correct vdjSQ_IB_correct vecRS vecIB vecIS vecRB xij_for_investorL plot_figuresL homedealerL etaL quantityrobusL dayL kappaL keepallL largetradeL tradesizeL zerobenefit fejL loyaltyL

%correct
vdjSQ_IS=vdjSQ_IS_correct;
vdjSQ_IB=vdjSQ_IB_correct;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%CREATE DATA VARIABLES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xij_for_investor homedealer eta quantityrobus day kappa Dmax Imax T keepall largetrade fej tradesize loyalty 
 
%Assing correct values (overwrite globals in the loaded files if necessary)
Imax                =9;
Dmax                =[1, 2, 3, 4, 5, 6, 7, 8, 9];
xij_for_investor    =xij_for_investorL;
homedealer          =homedealerL; 
eta                 =etaL;
quantityrobus       =quantityrobusL; 
day                 =dayL;
kappa               =kappaL;
keepall             =keepallL;
largetrade          =largetradeL;
tradesize           = tradesizeL;
fej                 = fejL;
loyalty             =loyaltyL;

Y=[YIB;YIS;YRB;YRS];
gvec = 1:4;

%sort on date, investor group
Y=sortrows(Y,[3,1]);
   
%Address for new dates (pool all groups)
new_period = zeros(length(unique(Y(:,3))),2);
count=2;
new_period(1, :)=[Y(1,3),1];
for i=2:size(Y,1)
    if Y(i,3)~=Y(i-1,3)
    new_period(count,:) = [Y(i,3), i];  
    count=count+1;
    end
end 
new_period_end = new_period(2:end,2)-1;
new_period = [new_period, [new_period_end;size(Y,1)]];

%Address for new dates (pool buy vs. sell)
new_period_side = zeros(length(unique(Y(:,3))),2+1);
count=2;
new_period_side(1, :)=[Y(1,3),1,Y(1,5)];
for i=2:size(Y,1)
    if Y(i,3)~=Y(i-1,3) || Y(i,5)~=Y(i-1,5) %same date, same side
    new_period_side(count,:) = [Y(i,3), i, Y(i,5)];
    count=count+1;
    end
end 

new_period_side_end = new_period_side(2:end,2)-1;
new_period_side = [new_period_side, [new_period_side_end;size(Y,1)]];
new_period_side = new_period_side(:,[1,2,4,3]);

T=size(new_period/2,1);
k=23+7;


%fill in c_t and value of the dealer for retail investors 
for t=1:T

    %Buyers
    ct  = unique(Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==1,k));
    if length(ct)>1 && sum(isnan(ct))==0
        disp('ct not unique')
    elseif sum(isnan(ct))>0
       disp('ct missing (was not estimated due to too few OTC trades)')
    end 
    
    if ~isnan(sum(unique(Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==1,k))))
        Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==2,k) = ct;
    else 
        Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==2,k) = NaN;
    end 

    %Sellers
    ct  = unique(Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==3,k));
    if length(ct)>1 && sum(isnan(ct))==0
        disp('ct not unique')
    elseif sum(isnan(ct))>0
       disp('ct missing (was not estimated due to too few OTC trades)')
    end 
    if ~isnan(sum(unique(Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==3,k))))
        Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==4,k) = ct;
    else 
        Y(Y(:,3)==Y(new_period(t,2),3) & Y(:,1)==4,k) = NaN;
    end   
end 
clear  ct


%Address for new dates (per group)
for g=gvec %group 1 IB, group 2 RB, group 3 IS, group 4 RS

    count=2;
    Y_g{g} = Y(Y(:,1)==g,2:end);
    new_period_g{g}(1, :)=[Y_g{g}(1,2),1];
    for i=2:size(Y_g{g},1)
        if Y_g{g}(i,2)~=Y_g{g}(i-1,2) 
        new_period_g{g}(count,:) = [Y_g{g}(i,2), i];  
        count=count+1;
        end
    end 
    
new_period_end_g{g} = new_period_g{g}(2:end,2)-1;
new_period_g{g}   = [new_period_g{g}, [new_period_end_g{g};size(Y_g{g},1)]];

T_g{g}=size(new_period_g{g},1);
end 

%Set option for maximization
options = optimoptions('fmincon' ,  'FiniteDifferenceType','central','Algorithm', 'interior-point', 'OptimalityTolerance',1e-40, 'StepTolerance', 1e-40, 'FunctionTolerance', 1e-40,   'ConstraintTolerance', 1e-40, 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6);  
lb =[0 0 0 0 0 0 0 0 0];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%STATUS QUO 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Status quo
SQ=1;                %means use the estimated value for the dealer
CF=0;                %means use esitmated cost of entry

for g=gvec

Yg          = Y_g{g};
new_periodg = new_period_g{g};
side        = unique(Yg(:,4));
retailer    = unique(Yg(:,3));

%to check
Yg_org    = Yg;
new_periodg_org = new_periodg;
if g==1 || g==2
    Yg          = Y_g{1};
    new_periodg= new_period_g{1};
else 
    Yg          = Y_g{3};
    new_periodg= new_period_g{3};
end 

%Create a Y-matrix with information for retailers other than qualitites of dealers 
%Note: in status quo fun the only thing that matters for retailers in the end is: Mhat, theta, dealer qualities and loyalty benefit
for t=1:size(new_periodg,1)
    idx_I      = new_periodg(t,2):new_periodg(t,3);
    idx_R    = new_periodg_org(t,2);
    Rep       = size(new_periodg(t,2):new_periodg(t,3),2);
    cols       = [3,6 8,9,10,11,12,13, 27,28,29] ;  
    Yg(idx_I, cols) = repmat(Yg_org(idx_R, cols),  Rep,1);
end 

%Value of dealer
if g==1 || g==2, vdjSQ = vdjSQ_IB;
elseif g==3 || g==4, vdjSQ = vdjSQ_IS;
end 

if retailer==0
    access=1;
elseif retailer==1
    access=0;
end 

[psiSQ,rhojSQ,sjSQ, sHjSQ, sNHjSQ, s0jSQ,s0HjSQ, s0NHjSQ, etaEjSQ,etaDjSQ,piDprimejSQ,dealer_idsSQ,uBSQ,uESQ,piBSQ,piESQ,YBSQ,YESQ,EcostSQ, ...
    q_data,rhoj_data, s0j_data,xij_data, Yno, EUno, EU2_SQ, psiSQ_opt,rhojSQ_opt, vdjSQout , WSQ, WxijSQ, WvdjSQ, WxijHSQ , WshrinkSQ, WSQwH, WxijSQwH, WvdjSQwH, WxijHSQwH , WshrinkSQwH] = status_quo_fun(SQ, CF, side, access, Yg, new_periodg, new_period, vdjSQ, 0, []);

%total yield, expected surplus and profit
YieldSQ     = YBSQ + YESQ;
EUSQ        = uBSQ + uESQ+ EcostSQ;
PiSQ        = piBSQ + piESQ; 

%sort according to dealers:
xij = xij_data;
med_xij=nanmedian(xij_data);
cols = [3,9,4,1,2,5,8,6,7];

if zerobenefit==0
    temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty)]]; 
    save(fullfile([foldermats,temp,'.mat']))
else 
    temp=['/counterfac_SQ_',[num2str(side),'_' num2str(retailer), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize), '_', num2str(fej), '_', num2str(loyalty), '_zerobenefit']]; 
    save(fullfile([foldermats,temp,'.mat']))
end   
end 
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%COUNTERFACTUALS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

kappavec =[0,1,0.1]; %share of retailers

for side=[1 -1]

%With flexible quotes
for CF=[-3 -2 -1 3 4 5 8 9 12 13 16 17 18]
    
    if CF==3|| CF==4 || CF==5 || CF==6 ||CF==16
        flexible=1;
    else 
        flexible=0;
    end 

    %Counterfactuals

    %CF= -3 no access to institutional investors  -- institutional investors only 

    %CF = -2 status quo but with q=v and sigma=0 -- institutional investors only 
    %CF = -1 status quo but with q=v  -- institutional investors only 
 
    %CF=0: at estimated costs, separate quotes and kappa=1 for R
    %CF=1: at estimate costs, separate quotes  and kappa=0 for I
    %CF=2: at 0 costs, separate quotes and  kappa=0 for I
    
    %CF=3: at estimated costs with the same quote (kappa=0.1 for R, (1-kappa)=0.9 for I)
    %CF=4: at 0 costs costs with the same quote (kappa=0.1 for R, (1-kappa)=0.9 for I)
    %CF=5: at at estimate costs -1.1 with the same quote (kappa=0.1 for R, (1-kappa)=0.9 for I)
    %CF=6: at 0 cost and 0 loyalty benefit  with the same quote (kappa=0.1 for R, (1-kappa)=0.9 for I)
           
    %CF=7: focus on institutional investors only. Shut down platform access at estimated costs (counterfac quotes can't be computed if no one is on the platform)
           
    %CF=8-11: Make the platform perfectly competitive (q=v and sigma = 0)
           %CF=8 at estimated cost, CF=9 eliminate fee, CF=10 eliminate cost and loyalty benefit, CF=11 eliminate all cost
    
    %CF=12-15: Assume dealers are truthful (q=v) but fix platform competititon (sigma as estimated)
           %CF=12 at estimated cost, CF=13 eliminate fee, CF=14 eliminate cost and loyalty benefit, CF=15 eliminate all cost
           
    %CF = 16-18: All on the platform
           %16 with quotes
           %17 with q=v but sigma>0
           %18 with q=v and sigma=0
                              
if CF~=7
    access=1;
elseif CF==7
    access=0;
end

if CF==0, kappa=1; end %This is used for the model fit table for platform access only
if CF<0 || CF==3 || CF==4|| CF==5 || CF==6 || CF>=8, kappa = 0.1;
elseif CF==1 || CF==2, kappa=0; %only institutional investors
elseif CF==7, kappa=0; %only institutional investors
end


%A) CF with access to retail investors
if CF>=0 
    
    %value of dealer
    if side==1,      vdjSQ = vdjSQ_IB;
    elseif side==-1, vdjSQ = vdjSQ_IS;
    end 

    [psicf21, rhojcf21, sjcf21, sHjcf21, sNHjcf21, s0jcf21, s0Hjcf21, s0NHjcf21, dealer_idscf21,qcf21,etaDjcf21,etaEjcf21, ...
    uBcf21,uEcf21,piBcf21,piEcf21,YBcf21,YEcf21,Ecostcf21, uE2cf21, ...
    Wcf21, Wxijcf21, Wvdjcf21, WxijHcf21, Wcf21_shrink,Wcf21wH, Wxijcf21wH, Wvdjcf21wH, WxijHcf21wH, Wcf21_shrinkwH, ...
    psicf21_g,rhojcf21_g,sjcf21_g,  sHjcf21_g, sNHjcf21_g, s0jcf21_g,  s0Hjcf21_g, s0NHjcf21_g, uBcf21_g,uEcf21_g,YBcf21_g, ...
    YEcf21_g,Ecostcf21_g, piBcf21_g,piEcf21_g, uE2cf21_g, ...
    psicf21_opt, rhojcf21_opt, psicf21_opt_g, rhojcf21_opt_g, ...
    Wcf21_g, Wxijcf21_g, Wvdjcf21_g, WxijHcf21_g, Wcf21_shrink_g, Wcf21_gwH, Wxijcf21_gwH, Wvdjcf21_gwH, WxijHcf21_gwH, Wcf21_shrink_gwH] ...
    = counterfac_per_side_fun(CF, side, access, Y,new_period,  new_period_side, vdjSQ, kappa, flexible);

    %total yield, expected surplus and profit
    Yield21     = YBcf21 + YEcf21; 
    EU21        = uBcf21 + uEcf21+ Ecostcf21; 
    Pi21          = piBcf21 + piEcf21; 


%B) CF with no access to retail investors but changes for institutionals 
else 

     retailer=0;
    if side==1
        Yg          = Y_g{1};
        new_periodg = new_period_g{1};
        vdjSQ = vdjSQ_IB;
    else 
        Yg          = Y_g{3};
        new_periodg = new_period_g{3};
        vdjSQ = vdjSQ_IS;
    end 
   
    if CF==-1 ||  CF==-2 %Access to institutional investors only, but make platform better 
        [psiSQ,rhojSQ,sjSQ, sHjSQ, sNHjSQ, s0jSQ,s0HjSQ, s0NHjSQ, etaEjSQ,etaDjSQ,piDprimejSQ,dealer_idsSQ,uBSQ,uESQ,piBSQ,piESQ,YBSQ,YESQ,EcostSQ, ...
        q_data,rhoj_data, s0j_data,xij_data, Yno, EUno, EU2_SQ, psiSQ_opt,rhojSQ_opt, vdjSQout , WSQ, WxijSQ, WvdjSQ, WxijHSQ , WshrinkSQ, WSQwH, WxijSQwH, WvdjSQwH, WxijHSQwH , WshrinkSQwH] = status_quo_fun(SQ, CF, side, 1, Yg, new_periodg, new_period, vdjSQ, 0, []);

    else %No access to institutional investors 
         [psiSQ,rhojSQ,sjSQ, sHjSQ, sNHjSQ, s0jSQ,s0HjSQ, s0NHjSQ, etaEjSQ,etaDjSQ,piDprimejSQ,dealer_idsSQ,uBSQ,uESQ,piBSQ,piESQ,YBSQ,YESQ,EcostSQ, ...
        q_data,rhoj_data, s0j_data,xij_data, Yno, EUno, EU2_SQ, psiSQ_opt,rhojSQ_opt, vdjSQout , WSQ, WxijSQ, WvdjSQ, WxijHSQ , WshrinkSQ, WSQwH, WxijSQwH, WvdjSQwH, WxijHSQwH , WshrinkSQwH] = status_quo_fun(SQ, CF, side, 0, Yg, new_periodg, new_period, vdjSQ, 0, []);
    end 

    YieldSQ     = YBSQ + YESQ;
    EUSQ        = uBSQ + uESQ+ EcostSQ;
    PiSQ        = piBSQ + piESQ;  

end 

%save to plot different investor groups in the same graph
if zerobenefit==0
    temp=['/counterfac_',[num2str(side), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize),'_' num2str(CF),'_',num2str(flexible), '_', num2str(fej), '_', num2str(loyalty)]];
    save(fullfile([foldermats,temp,'.mat']))
else 
    temp=['/counterfac_',[num2str(side), '_',num2str(day), '_', num2str(homedealer), '_', num2str(keepall),'_', num2str(quantityrobus),'_', num2str(largetrade), '_', num2str(tradesize),'_' num2str(CF),'_',num2str(flexible), '_', num2str(fej), '_', num2str(loyalty),  '_zerobenefit']];
    save(fullfile([foldermats,temp,'.mat']))
end 

end 
end 

 




